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Gas hydrates are systems of prime importance. In particular, hydrogen hydrates are potential materials of 
icy satellites and comets, and maybe used for hydrogen storage. We explore the H 2 0-H 2 system at pressures 
in the range 0-100 GPa with ab initio variable-composition evolutionary simulations. According to our 
calculation and previous experiments, the H 2 0-H 2 system undergoes a series of transformations with 
pressure, and adopts the known open-network clathrate structures (sll, C 0 ), dense "filled ice" structures (Ci, 
C 2 ) and two novel hydrate phases. One of these is based on the hexagonal ice framework and has the same 
H 2 0 : H 2 ratio (2 : 1) as the C 0 phase at low pressures and similar enthalpy (we name this phase Ih-C 0 ). The 
other newly predicted hydrate phase has a 1 : 2 H 2 0 : H 2 ratio and structure based on cubic ice. This phase 
(which we name C 3 ) is predicted to be thermodynamically stable above 38 GPa when including van der 
Waals interactions and zero-point vibrational energy, and explains previously mysterious experimental 
X-ray diffraction and Raman measurements. This is the hydrogen-richest hydrate and this phase has a 
remarkable gravimetric density (18 wt.%) of easily extractable hydrogen. 

Molecular compounds (cocrystals) of water ice (H 2 0) and hydrogen (H 2 ) are known to form clathrate 
structures with hydrogen molecules encapsulated as guests in the host sublattice formed by water 
molecules. Hydrogen hydrates, as environmentally clean and efficient hydrogen storage materials, have 
excited significant interest. Extensive literature exists from both experimental 1 " 12 and theoretical 1314 sides. Aside 
from the H 2 molecules, many other small molecules can be encapsulated as guests in clathrates at elevated 
pressure (e.g., noble gases, nitrogen, oxygen, methane, chlorine). (See Ref. [15] and references therein) 
Hydrogen hydrates are important as potentially major materials of icy satellites and comets, and hydrogen 
storage materials. 

Twenty years after the first report of the formation of two filled-ice hydrogen hydrates by Vos et al 1 , four 
hydrogen hydrate forms are known to exist at elevated pressures. Two of the hydrogen hydrates are clathrates, 
denoted as clathrate structure II (sll) 3 ' 5 and compound 0 (C 0 ) 12 , the other two are filled ice hydrates, compound 1 
(Ci) and compound 2 (C 2 ) 1,2 . The sll clathrate hydrate was synthesized under pressures of 180 to 220 MPa at 
300 K, and its structure was shown to contain 48 hydrogen molecules and 136 water molecules in the unit cell 3 . 
The C 0 clathrate was recently found to be stable near 0.5 GPa and to have the composition 2H 2 0 : 1H 2 and a 
trigonal quartz-like structure 12 . The water molecules in the C 0 structure are arranged in a totally new way, 
different from the known ices or ice sublattices in hydrates structures. This structure has space group P3 2 21, 
but this could possibly go as low as P3 2 , depending on how the hydrogens are arranged 12 . 

At higher pressures, clathrates give way to denser structures of the filled ice type. The Q and C 2 phases are 
formed at 0.36-0.9 GPa and —2.4 GPa, respectively 1 ' 211 . The C x hydrate has a water host framework based on ice- 

II and a 6 : 1 water to hydrogen ratio. C 2 has a 1 : 1 ratio of water to hydrogen and is composed of water molecules 
in the "cubic ice" (ice-Ic) framework and rotationally disordered hydrogen molecules 15 . Recent experiments 7-10 
indicate that the C 2 hydrate undergoes a structural transformation from cubic to tetragonal phase at around 10- 
20 GPa, with an increasing difference in the unit cell axes, and then transforms to another high-pressure phase 
near —45 GPa. This high-pressure phase is maintained up to at least 80 GPa but its structure is not fully resolved. 
Given the difficulties in characterization of the chemical composition and crystal structure of these hydrates, and 
believing that new phases are likely to exist, we decided to perform a computational search to revisit the H 2 0-H 2 
system under pressure. 
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Figure 1 | Convex hull diagram for H 2 0-H 2 system at selected pressures and zero temperature. This figure shows the enthalpy of formation (in eV/ 
molecule) of molecular compounds from pure H 2 0 ice and H 2 . The red and yellow circles represent the C 0 and Ih-C 0 phases, respectively. The green star 
represents the sll structure. 



Results 

Using the evolutionary algorithm USPEX 16 " 20 , we explored all pos- 
sible stable phases in the H 2 0-H 2 system (See Methods). 
Remarkably, we have found two novel filled ice hydrogen hydrates, 
and all known hydrogen hydrates (except the sll structure, because of 
the very large number of molecules in its unit cell). Thus, at pressures 
in the range 0-2 GPa, the sll structure was input separately in order 
to calculate stability ranges of phases in the H 2 0-H 2 system. Fig. 1 
shows the convex hull diagram for the H 2 0-H 2 system. 

With including the van der Waals (vdW) dispersion forces (See 
Methods), our results are generally in very good agreement with 
experiments, but with several novel aspects. At 0 GPa, the C 0 , C x 
and a novel hydrogen hydrate phase are found stable or nearly stable 
in the H 2 0-H 2 system, while the sll phase is metastable (—0.013 eV/ 
molecule less stable than the mixture of stable compounds C 0 and 
Ci). The structure of the novel hydrogen hydrate is based on the 
framework of hexagonal ice (ice-Ih), with two hydrogen molecules 
hosted inside channels running along the hexagonal axis (Fig. 2a). It 
has a 2 : 1 ratio of water to hydrogen, same as C 0 , and space group Cc. 
We name it Ih-C 0 to distinguish from C 0 . The enthalpy of the Ih-C 0 
phase is close to C 0 , and is slightly lower at pressures above —0.4 GPa 
(Fig. S2 in Supplementary Information). At 1.5 GPa, in addition to 
the C 0 , Ih-C 0 and C x phases, the hydrate phase C 2 with an ice-Ic 
framework structure becomes stable. 

At pressures above 2 GPa, the C 0 and Ih-C 0 phases are calculated 
to be above the convex hull, indicating that these phases become 
unstable against decomposition into C x and C 2 . Above 3.5 GPa, 
the Ci phase will also become unstable, and the C 2 phase will remain 
the only stable hydrate. For hydrate C 2 , USPEX calculations uncov- 
ered at least four typical energetically favorable candidate structures 
at different pressures, I^/amd, J>na2 x and l\ x md (Fig. S3, in 

Supplementary Information). These structures differ in the orienta- 




Figure 2 | Crystal structures of Ih-C 0 and C 3 phases, cage structures in C 2 
and C 3 phases, (a) Ih-C 0 structure at 0.5 GPa, (b) C 3 structure at 30 GPa, 
(c) cages formed by water molecules in C 3 at 100 GPa, the hydrogen 
molecules are located at the center of each chair-like H-O ring, (d) cages in 
"filled ice-Ic" C 2 , hydrogen molecules are in the center of the cage. Large 
red and small blue spheres are O and H atoms in water molecules, 
respectively; the yellow spheres represent the H 2 molecules in (a) and (b), 
and represent H atoms in (c) and (d). Red dashed lines represent hydrogen 
bonds. 
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Figure 3 | Phase diagram of the H 2 0-H 2 system. The stability ranges of C 2 and C 3 phases are calculated with and without ZPE effect. The solid 
orange line represents extra stability range added due to ZPE, the dashed orange line represents regions that become unstable after inclusion of the 
ZPE. 



tions of water and hydrogen molecules, which is similar to Ref. [15]. 
The C 2 phase will lose stability at —14 GPa, which is much lower 
than 40 GPa suggested in the previous study 710 . We explain this by 
metastable persistence of C 2 up to the pressure of 40 GPa. Between 
14-28 GPa, there are, unexpectedly, no thermodynamically stable 
hydrates. 

Near 30 GPa, another novel H 2 0-H 2 phase is found to be stable at 
zero temperature. It has a 1 : 2 water to hydrogen ratio, and net 
composition H 6 0. This novel hydrogen hydrate, which we name 
C 3 , has the highest hydrogen concentration among all hydrogen 
hydrates. If it can be synthesized at low pressures, it would be an 
attractive hydrogen storage material, having 18 wt.% concentration 
of easily separable (non-water) hydrogen. The C 3 structure has space 
group P4 X and is also based on the framework of ice-Ic (Fig. 2b), 
similar to low-pressure hydrate C 2 . The unit cell of C 3 contains four 
water molecules, the H 2 molecules are located at the center of chair- 
like H-O rings (formed by six oxygen and six hydrogen atoms) that 
form faces of the cage, as shown in Fig. 2c. Differently, in the C 2 
hydrate, the H 2 molecules are in the center of the water cages 
(Fig. 2d). According to our calcualtions, the C 3 phase will remain 
stable up to at least 120 GPa. 

Our theoretical calculations indicate that the H 2 0-H 2 system 
contains several stable phases, including open-network clathrate 
structures (C 0 ) and dense filled ice phases (Ih-C 0 , C l5 C 2 and C 3 ). 
The C 0 phase is predicted to be stable at pressures below 1.5 GPa, 
which is close to the experimental result (below 0.8 GPa 12 ). The Ci 
phase is predicted to be stable at pressures below 3.5 GPa, also close 
to the experimentally determined transition pressure of 2.5 GPa 1 . 
The zero-point vibration energy (ZPE) significantly affects the rela- 
tive stability of hydrogen- rich structures 21 . We have estimated the 
ZPE within the quasi-harmonic approximation 22 to refine the 
stability ranges of C 2 and C 3 phases above 10 GPa. When consider- 
ing the ZPE, the stability field of the C 2 phase expands up to 
— 19 GPa, but this phase remains dynamically stable, and thus can 
exist as a metastable material at pressures of at least 60 GPa (Fig. S9 
and Fig. S10 in Supplementary Information). 

The C 3 phase starts to be energetically favorable above —38 GPa 
when including ZPE, as shown in Fig. 3. Thus, the novel C 3 phase can 



be synthesized in hydrogen -rich conditions at pressures starting 
from 38 GPa. This theoretical value agrees well with the experiment- 
ally observed formation at 45-50 GPa of a mysterious phase of 
unknown composition 710 . As shown in Fig. 4, calculations of the 
Raman shift 23 reveal differences between the C 2 and C 3 phases in 
the H 2 -D 2 0 system. The Raman shift of C 3 phase, rather than an 
amorphous phase, agrees very well with the lower Raman frequencies 
of the vibron for the hydrogen molecules observed in Ref. [9] . The 
black rhombi in Fig. 4 indicate that some of the H 2 -D 2 0 C 3 sample 
encountered decomposition when quenched to low pressure. The 
variation of lattice parameters of the ice host structure in hydrates 
with pressure, revealed by our theoretical calculations, also agrees 
well with experimental XRD results 10 . At 55 GPa, our calculation 
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Figure 4 | Variations of the Raman shift of the vibron for the H 2 
molecules with pressure from experimental data in Ref. [9] and our 
theoretical calculations. The red and black symbols are the experimental 
data for H 2 vibrons in the H 2 -D 2 0 sample. The blue open circles and 
squares indicate the Raman shift calculation for C 2 and C 3 phases of H 2 - 
D 2 0 system, respectively. 
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Figure 5 | Internal energy of the C 2 and C 3 phases relative to the 
isochoric mixture of H 2 0 and H 2 . Green lines represent the energy of the 
hydrate phases; red lines represent the energy of the isochoric mixture of 
ice-VIII and H 2 -I phases. 

gives lattice parameter of C 3 phase a = b = 4.00 A and c = 5.67 A, 
corresponding to cubic ice sublattice with periodicity 5.67 A, 
whereas experiment gives —5.5 A 10 . 

At low pressure, the C 2 adopts a "cubic ice" host structure and 
then transforms to a "tetragonal" one at —20 GPa 10 (Fig. S4 in 
Supplementary Information). When forming the C 3 phase at 
increased pressure and in excess of H 2 , the ice host structure trans- 
forms to the "cubic ice" again. The change from tetragonal to "cubic" 
structure occurs before H-bond symmetrization transition happens 
in "tetragonal" type C 2 at —55 GPa. Thus, such structural trans- 
formation is unrelated to symmetrization of the H -bonds, but comes 
from the emergence of the C 3 phase. For the hydrate C 3 , the H-bond 
symmetrization is predicted to occur at —120 GPa (Fig. S5 in 
Supplementary Information), which is close the theoretical H-bond 
symmetrization pressure in ice-VII 14 . 

Discussion 

The C 2 and C 3 hydrates have a similar ice host framework, but the 
different numbers of hydrogen molecules and their different loca- 
tions and orientations bring huge differences in phase stability range. 
In the C 2 phase, hydrogen molecules stay in the centers of cages 
formed by water molecules in contrast to C 3 phase, where they are 
located at the faces of the cages. To clarify the causes of stability of 
hydrogen hydrates, we used Bader analysis 24,25 , and focused on the C 2 
and C 3 phases (Fig. S6 in Supplementary Information). We found a 
very small charge transferred from H 2 to water molecules, so that the 
H 2 molecules are slightly positively charged, and H 2 0 molecules 
carry a slight negative charge. The magnitudes of these charges are 
— 10" 3 -10" 2 per molecule. This suggests that interactions between 
these molecules are almost purely steric, mainly related to packing 
density and shapes of the molecules. Comparing Bader volumes of 
the H 2 0 and H 2 molecules in the hydrates and in pure H 2 0 and H 2 , 
we see that water molecules occupy slightly larger volume in the 
hydrates, whereas hydrogen molecules occupy much less space in 
C 3 hydrate than in pure H 2 - this leads to net densification, stabil- 
izing this phase in a wide pressure range. For the C 2 hydrate, the H 2 
molecules have lower volume than in pure H 2 only at pressures below 
— 10 GPa, which explains its instability at higher pressures. 

Having considered the Py-term in the enthalpy (H = E + PV), to 
get additional insight, we turned to the internal energy E and its 
changes when the H 2 and H 2 0 molecules are placed from the hydrate 



into pure H 2 and H 2 0 phases, while keeping molecular volumes fixed 
at their values in the hydrate (Fig. 5). This energy characterizes the 
net balance of the vdW attraction and steric repulsion between the 
molecules: this net effect is very small in the C 3 phase (slightly desta- 
bilizing below —30 GPa and slightly stabilizing above —30 GPa). 
The remarkably wide stability field of the C 3 phase is therefore mostly 
due to its high density and only to a small extent to more favorable 
intermolecular interactions. A much more interesting picture is 
observed for the C 2 phase (Fig. 5a): we find its slight energetic sta- 
bilization below —15 GPa, and an increasingly large destabilization 
at higher pressures. This explains why C 2 is unstable at high pres- 
sures, and furthermore, it is clear that the increasing energetic instab- 
ility of the C 2 phase is responsible for the displacive phase transition, 
metastably occurring on overcompression and transforming the 
cubic H 2 0 host sublattice into tetragonal, to enable better packing 
of the molecules. 

Our calculations found that a C 3 -type phase is stable in the H 2 0- 
He system at 8-75 GPa (without including zero-point energy), and 
this phase is denser than the mixture of H 2 0 and He. On the other 
hand, no such phase was found in the H 2 0-Ne system, and indeed 
the C 3 phase is not packing- efficient in this system (Fig. S7 and Fig. 
S8 in Supplementary Information). He and Ne are both chemically 
extremely inert, their almost only differences are size and (here 
insignificant) mass. Stability of He-C 3 and instability of Ne-C 3 
hydrates reinforce our conclusion made for the H 2 0-H 2 system, that 
stability of this novel phase comes not from specific bonding inter- 
actions between the molecules, and not even their shapes, but is 
mostly due to their very efficient packing. 

In summary, using the evolutionary algorithm USPEX, we 
explored the H 2 0-H 2 system at pressures of up to 100 GPa. 
Stoichiometries and stability fields of H 2 0-H 2 hydrate phases have 
been studied. A series of pressure-induced transformations found by 
theory closely coincides with experimental data, but also new insight 
was obtained. A novel Ih-C 0 structure is predicted to have a very 
close enthalpy to the recently discovered C 0 structure. At pressures 
above 38 GPa, novel hydrogen hydrate C 3 , based on cubic ice Ic, is 
predicted to be stable. With stoichiometry H 2 0:2H 2 , this is the 
hydrogen -richest hydrate known to date. With gravimetric density 
of easily removable hydrogen (18 wt.%), this is a promising hydro- 
gen storage material that can find practical applications if its syn- 
thesis pressure can be decreased. 

Methods 

Crystal structure prediction. Predictions were done using the USPEX code in its 
variable-composition mode at several pressures (0, 1, 2, 5, 10, 20, 50 and 100 GPa) 
and zero temperature. A number of studies illustrate the power of the USPEX 
method 26 " 28 . We have done two types of variable-composition structure predictions in 
searching for all stable phases in the H-O system: (1) In the H-O system, assembling 
the structures from atoms, and (2) In the H 2 0-H 2 system with giving H 2 0 and H 2 
molecules as structure building blocks. We have found that, even in the H-O system, 
all low-enthalpy states at pressures of our interest (<120 GPa) are actually made of 
well-defined H 2 0 and H 2 molecules. This allowed us to focus on molecular-type 
calculations, capable of efficiently dealing with large systems, without loss of rigor. 

Given molecular nature of all stable and nearly stable compounds in this system, we 
searched for the packing of well-defined H 2 0 and H 2 molecules (rather than H and O 
atoms), by applying the specially designed constrained global optimization algorithm, 
considering structures with up to 24 molecules (i.e. up to 72 atoms) per primitive unit 
cell. 

DFT calculations. Structure relaxations were done using density functional theory 
(DFT) within van der Waals (vdW) functional optB88-vdW 29 in the framework of the 
all-electron projector augmented wave (PAW) 30 method as implemented in the 
VASP 31 code. The plane wave kinetic energy cutoff of 600 eV and Gamma-centered 
&>point meshes with the reciprocal space resolution of 2n X 0.05 A were used. Having 
identified the most stable compositions and candidate structures, we relaxed them at 
pressures from 1 atm to 120 GPa with an even higher cutoff of 800 eV to refine their 
thermodynamic properties and stability fields. Structure relaxations proceeded until 
net forces on atoms were below 1 meV/A, which gave us enthalpies converged to 
better than 1 meV/atom. 

It is expected that the relative contribution of hydrogen bonding (H-bonding) and 
van der Waals (vdW) dispersion forces has a significant impact on the phase trans- 
ition pressures and cohesive properties of the various crystalline ice phases 32 . This is 
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also confirmed by our calculations for the phase transition pressures of ice phases 
from optB-88-vdW, GGA 33 calculations and experiments (Fig. SI in Supplemental 
Information). Thus, all calculations included the vdW functional to treat the vdW 
forces, unless stated otherwise. 
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